//
//
//---------------------  2D  ---------------------
//
//       Staggered Mesh for u-vel and v-vel
//
//       0       1       2       3       4          
//
//5      >       >       >       >       > 
//       |       |       |       |       |              
//   ^---+---^---+---^---+---^---+---^---+---^  4       Mesh for scalar fields
//       |       |       |       |       | 
//4      >   o   >   o   >   o   >   o   >               5  x-x-+-x-+-x-+-x-x
//       |       |       |       |       |               4  x o | o | o | o x
//   ^---+---^---+---^---+---^---+---^---+---^  3           +---+---+---+---+
//       |       |       |       |       |               3  x o | o | o | o x
//3      >   o   >   o   >   o   >   o   >                  +---+---+---+---+
//       |       |       |       |       |               2  x o | o | o | o x
//   ^---+---^---+---^---+---^---+---^---+---^  2           +---+---+---+---+
//       |       |       |       |       |               1  x o | o | o | o x
//2      >   o   >   o   >   o   >   o   >               0  x-x-+-x-+-x-+-x-x
//       |       |       |       |       |                  
//   ^---+---^---+---^---+---^---+---^---+---^  1           0 1   2   3   4 5 
//       |       |       |       |       |                 
//1      >   o   >   o   >   o   >   o   >                 o central node
//       |       |       |       |       |                 x boundary node
//   ^---+---^---+---^---+---^---+---^---+---^  0          > u velocity 
//       |       |       |       |       |                 ^ v velocity
//0      >       >       >       >       >    
//   0       1       2       3       4       5       
//                                                  
//                      v(i,j) 
//                  |     n     |               
//                --|-----^-----|--             
//                  |           |               
//                  |           |               
//     u(i-1,j) = w >     o     > e = u(i,j)
//                  |   (i,j)   |
//                  |           |
//                --|-----^-----|--
//                  |     s     |
//                     v(i,j-1)
//     


namespace Tuna {

template<class Tprec, int Dim>
inline bool Quick_CoDi<Tprec, Dim>::calcCoefficients1D() 
{
    prec_t Gamma_dx = Gamma / dx;
    prec_t dx_dt = dx / dt;
    prec_t ce, cem, cep, cw, cwm, cwp;
    aE = 0.0; aW = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int i =  bi; i <= ei; ++i) {
	ce = u(i  );
	cw = u(i-1);

	// OJO : los terminos cem y cwm son POSITIVIOS!!!
	// Checar los casos 2D y 3D!!!
	    
	if ( ce > 0 ) { cem = 0.0; cep = ce * 0.125; }
	else {          cem = ce * 0.125; cep = 0.0; }
	
	if ( cw > 0 ) { cwm = 0.0; cwp = cw * 0.125; }
	else {          cwm = cw * 0.125; cwp = 0.0; }

	aE (i) = Gamma_dx - ce * 0.5 + cep - 2 * cem - cwm;
	aW (i) = Gamma_dx + cw * 0.5 + 2 * cwp - cwm + cep; 
	aP (i) = aE (i) + aW (i) + dx_dt + cem - cwp;
		+ (ce - cw) ;
// Term (ce - cw) represents the discretizated continuity equation, and
// must be equal to zero when that equation is valid, so I can avoid to
// calculate this term for efficiency.
	    
	sp (i) = phi_0(i) * dx_dt ;

	if (i <= ei-1)        { sp (i) += cem * phi_0(i+2); }
	else if ( i == ei) { sp (i) += cem * phi_0(i+1); }
	
	if (i >= bi+1)        { sp (i) -= cwp * phi_0(i-2); }
	else if ( i == bi) { sp (i) -= cwp * phi_0(i-1); }       	
    }

    /*** *
    // Dirichlet Boundary Condition on LEFT WALL
    prec_t phi_A = 1;
    aP(bi) += aW(bi) + Gamma_dx;
    aE(bi) += Gamma_dx / 3.0;  
    sp(bi) += (2 * aW(bi) + 2.0 * Gamma_dx / 3.0) * phi_A;
    aW(bi) = 0.0;

    /*** /
    // Dirichlet Boundary Condition on LEFT WALL
    prec_t phi_B = 0;
    aP(ei) += aE(ei) + Gamma_dx;
    aW(ei) += Gamma_dx / 3.0;  
    sp(ei) += (2 * aE(ei) + 2.0 * Gamma_dx / 3.0) * phi_B;
    aE(ei) = 0.0;

    /* ****/

    applyBoundaryConditions1D();
    return 0;
}
          
template<class Tprec, int Dim>
inline bool Quick_CoDi<Tprec, Dim>::calcCoefficients2D() 
{
    prec_t Gdy_dx = Gamma * dy / dx;
    prec_t Gdx_dy = Gamma * dx / dy;
    prec_t dxy_dt = dx * dy / dt;
    prec_t ce, cem, cep, cw, cwm, cwp;
    prec_t cn, cnm, cnp, cs, csm, csp;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int i =  bi; i <= ei; ++i)
	for (int j = bj; j <= ej; ++j)
	{
	    ce = u(i, j) * dy;
	    cw = u(i-1, j) * dy;
	    cn = v(i, j) * dx;
	    cs = v(i, j-1) * dx;
	    
	    if ( ce > 0 ) { cem = 0.0; cep = ce * 0.125; }
	    else {          cem = ce * 0.125; cep = 0.0; }
	    
	    if ( cw > 0 ) { cwm = 0.0; cwp = cw * 0.125; }
	    else {          cwm = cw * 0.125; cwp = 0.0; }

	    if ( cn > 0 ) { cnm = 0.0; cnp = cn * 0.125; }
	    else {          cnm = cn * 0.125; cnp = 0.0; }
		
	    if ( cs > 0 ) { csm = 0.0; csp = cs * 0.125; }
	    else {          csm = cs * 0.125; csp = 0.0; }

	    aE (i,j) = Gdy_dx - ce * 0.5 + cep - 2 * cem - cwm;
	    aW (i,j) = Gdy_dx + cw * 0.5 + 2 * cwp - cwm + cep; 
	    aN (i,j) = Gdx_dy - cn * 0.5 + cnp - 2 * cnm - csm;
	    aS (i,j) = Gdx_dy + cs * 0.5 + 2 * csp - csm + cnp;
	    aP (i,j) = aE (i,j) + aW (i,j) + aN (i,j) + aS (i,j) + dxy_dt +
		cem - cwp + cnm - csp;
//		+ (ce - cw) + (cn - cs) ;
// Term (ce - cw) + (cn - cs) represents the discretizated continuity equation, and
// must be equal to zero when that equation is valid, so I can avoid
// this term for efficiency.
	    
	    sp (i,j) = phi_0(i,j) * dxy_dt ;

	    if (i <= ei-1)        { sp (i,j) += cem * phi_0(i+2,j); }
	    else if ( i == ei) { sp (i,j) += cem * phi_0(i+1,j); }

	    if (i >= bi+1)        { sp (i,j) -= cwp * phi_0(i-2,j); }
	    else if ( i == bi) { sp (i,j) -= cwp * phi_0(i-1,j); }

	    if (j <= ej-1)        { sp (i,j) += cnm * phi_0(i,j+2); }
	    else if ( j == ej) { sp (i,j) += cnm * phi_0(i,j+1); }

	    if (j >= bj+1)        { sp (i,j) -= csp * phi_0(i,j-2); }
	    else if ( j == bj) { sp (i,j) -= csp * phi_0(i,j-1); }

	}    
    applyBoundaryConditions2D();
    return 0;
}


//
//---------------------  3D  ---------------------
//
template<class Tprec, int Dim>
inline bool Quick_CoDi<Tprec, Dim>::calcCoefficients3D()
{
    prec_t dyz = dy * dz, dyz_dx = Gamma * dyz / dx;
    prec_t dxz = dx * dz, dxz_dy = Gamma * dxz / dy;
    prec_t dxy = dx * dy, dxy_dz = Gamma * dxy / dz;
    prec_t dxyz_dt = dx * dy * dz / dt;
    prec_t ce, cem, cep, cw, cwm, cwp;
    prec_t cn, cnm, cnp, cs, csm, csp;
    prec_t cf, cfm, cfp, cb, cbm, cbp;
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aF = 0.0; aB = 0.0; aP = 0.0; 
    sp = 0.0;

    for (int k = bk; k <= ek; ++k)
	for (int i =  bi; i <= ei; ++i)
	    for (int j = bj; j <= ej; ++j)
	    {
		ce = u(i  , j, k) * dyz;
		cw = u(i-1, j, k) * dyz;
		cn = v(i, j  , k) * dxz;
		cs = v(i, j-1, k) * dxz;
		cf = w(i, j, k  ) * dxy;
		cb = w(i, j, k-1) * dxy;

		if ( ce > 0 ) { cem = 0.0; cep = ce * 0.125; }
		else {          cem = ce * 0.125; cep = 0.0; }
		
		if ( cw > 0 ) { cwm = 0.0; cwp = cw * 0.125; }
		else {          cwm = cw * 0.125; cwp = 0.0; }

		if ( cn > 0 ) { cnm = 0.0; cnp = cn * 0.125; }
		else {          cnm = cn * 0.125; cnp = 0.0; }
		
		if ( cs > 0 ) { csm = 0.0; csp = cs * 0.125; }
		else {          csm = cs * 0.125; csp = 0.0; }

		if ( cf > 0 ) { cfm = 0.0; cfp = cf * 0.125; }
		else {          cfm = cf * 0.125; cfp = 0.0; }
		
		if ( cb > 0 ) { cbm = 0.0; cbp = cb * 0.125; }
		else {          cbm = cb * 0.125; cbp = 0.0; }

		aE (i,j,k) = dyz_dx - ce * 0.5 + cep - 2 * cem - cwm;
		aW (i,j,k) = dyz_dx + cw * 0.5 + 2 * cwp - cwm + cep; 
		aN (i,j,k) = dxz_dy - cn * 0.5 + cnp - 2 * cnm - csm;
		aS (i,j,k) = dxz_dy + cs * 0.5 + 2 * csp - csm + cnp;
		aF (i,j,k) = dxy_dz - cf * 0.5 + cfp - 2 * cfm - cbm;
		aB (i,j,k) = dxy_dz + cb * 0.5 + 2 * cbp - cbm + cfp; 
		aP (i,j,k) = aE (i,j,k) + aW (i,j,k) 
		    + aN (i,j,k) + aS (i,j,k)
		    + aF (i,j,k) + aB (i,j,k) + dxyz_dt
		    + cem - cwp + cnm - csp + cfm - cbp;
//		    + (ce - cw) + (cn - cs) + (cf - cb);
// Term (ce - cw) + (cn - cs) + (f -cb) represents the discretizated continuity 
// equation, and must be equal to zero when that equation is valid, so I can avoid
// this term for efficiency.

		sp (i,j,k) = phi_0(i,j,k) * dxyz_dt ;

		if (i <= ei-1)        { sp (i,j,k) += cem * phi_0(i+2,j,k); }
		else if ( i == ei) { sp (i,j,k) += cem * phi_0(i+1,j,k); }
		
		if (i >= bi+1)        { sp (i,j,k) -= cwp * phi_0(i-2,j,k); }
		else if ( i == bi) { sp (i,j,k) -= cwp * phi_0(i-1,j,k); }

		if (j <= ej-1)        { sp (i,j,k) += cnm * phi_0(i,j+2,k); }
		else if ( j == ej-1) { sp (i,j,k) += cnm * phi_0(i,j+1,k); }
		
		if (j >= bj+1)        { sp (i,j,k) -= csp * phi_0(i,j-2,k); }
		else if ( j == bj) { sp (i,j,k) -= csp * phi_0(i,j-1,k); }

		if (k <= ek-1)        { sp (i,j,k) += cfm * phi_0(i,j,k+2); }
		else if ( k == ek) { sp (i,j,k) += cfm * phi_0(i,j,k+1); }
		
		if (k >= bk+1)        { sp (i,j,k) -= cbp * phi_0(i,j,k-2); }
		else if ( k == bk) { sp (i,j,k) -= cbp * phi_0(i,j,k-1); }
	}    
    applyBoundaryConditions3D();
    return 0;
}

} // Tuna namespace
















